/*************************************************************************
         Gravitational Accelerations Inverter and Analyzer (GAIA)
**************************************************************************
Purpose: GAIA analyzes GRACE acceleration data by defining a grid on
         the surface of the Earth and solving for various parameters at
         those grid points, using the method specified in the function
         create_options(). For example, Newtonian gravity can be used
         to invert the acceleration data for mass concentrations (mascons)
         on the Earth's surface that can vary with time sinusoidally
         (to recover tides and annual hydrological cycles, etc) or
         vary as a polynomial of time (constant, trend, quadratic, etc).
**************************************************************************
Input:   Command line arguments can be used to add a description to the
         folder that this program will create and fill with data.
         Recommended use: ./cgaia.sh "your description here"
**************************************************************************
Output:  Output files are created for create_plots.pro which reads them
         and makes plots (a map, a 2D time series plot, etc).
         These plots are placed in a folder named using the date and the
         command line arguments (in GAIAFOLDER/output/maps/).
         "Most" source code is compressed and saved in that folder too.
**************************************************************************
Copyright 2010 Bryan Killett

This file is part of GAIA.

GAIA is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

GAIA is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with GAIA. If not, see <http://www.gnu.org/licenses/>.
*************************************************************************/
//Header file "definitions.hpp" declares many functions, constants, datatypes.
//In order to use the functions declared in definitions.hpp, must link
//with object code that includes functions.cpp.
#include "definitions.hpp"

void create_options(string &options_file, int argc, char *argv[]){
/**********************************************************************
Purpose: This function sets various options for GAIA using hard-coded
         values. The chosen options are written to disk and reloaded in
         load_options().
***********************************************************************
Input:   (types ending in "_s" are defined in definitions.hpp)
         output_file - filename of options.txt
         argc        - number of words in command line string
         argv        - those words, used to make output folder name
***********************************************************************
Output:  (types ending in "_s" are defined in definitions.hpp)
         writes text file options.txt to disk.
**********************************************************************/
input_s input;
ancillary_s ancillary;
grid_options_s region_options;
grid_s grid;
model_options_s model_options;
analysis_options_s analysis_options;
plot_options_s plot_options;
insitu_s insitu;//For tide gauge data, etc.

////////////////////////////////////////////////////////////////////////////
//Analysis options that are often changed.
////////////////////////////////////////////////////////////////////////////

//input.options.filename="input_7_yrs_20s_smoothing.txt";
//input.options.filename="input_deg150_2008_20s_smoothing_global.txt";
//input.options.filename="input_deg90_2008_20s_smoothing_50N.txt";
//input.options.filename="input_deg90_2007_20s_smoothing_50N.txt";
//input.options.filename="input_deg150_2008_20s_smoothing_50N.txt";
//input.options.filename="input_and_orbit_deg90_minus_deg150_2008_20s_smoothing.txt";
//input.options.filename="input_and_orbit_deg90_minus_deg150_OCT2008.txt";
//input.options.filename="input_and_orbit_deg150_2008_20s_smoothing_50N.txt";
//input.options.filename="input_and_orbit_7yrs_fixed_orbit.txt";
//input.options.filename="input_and_orbit_land_and_ocean_edited.txt";
//input.options.filename="input_and_orbit_land_and_ocean_edited_cecilia_times.txt";
//input.options.filename="input_and_orbit_land_and_ocean_last_500k.txt";
//input.options.filename="input_and_orbit_last_500k.txt";
//input.options.filename="input_and_orbit_noise_from_land_ocean.txt";
//input.options.filename="input_and_orbit_sine_50s_from_land_ocean.txt";
//input.options.filename="input_and_orbit_noise_plus_sine_50s.txt";
//input.options.filename="input_and_orbit_land_and_ocean_only_1_of_10_pts.txt";
//input.options.filename="input_and_orbit_land_and_ocean_ascending_only.txt";
//input.options.filename="input_and_orbit_land_and_ocean_descending_only.txt";
//input.options.filename="input_and_orbit_short2.txt";//150000 pts.
//input.options.filename="input_and_orbit_land_ocean_short2.txt";//150000 pts.
//input.options.filename="input_and_orbit_land_and_ocean_edited_1_of_5.txt";
input.options.filename="input_and_orbit_short3.txt";
//input.options.filename="input_and_orbit_75N.txt";
//input.options.filename="input_and_orbit_75N_250k_pts.txt";
//input.options.filename="input_and_orbit_sim_type_76_75N_250k_pts_hmc1.txt";
//input.options.filename="input_and_orbit_sim_type_76_75N_250k_pts_hmc2_nmax200.txt";
//input.options.filename="input_and_orbit_sim_type_76_75N_250k_pts_hmc202_supp3.txt";
//input.options.filename="input_and_orbit_sim_type_76_75N_250k_pts_hmc202_supp15.txt";
//input.options.filename="input_and_orbit_hudson_swaths.txt";
//input.options.filename="input_and_orbit_sim_type_200_supp3.txt";
//input.options.filename="input_and_orbit_sim_type_119_minus_122_and_125.txt";
//./cgaia.sh "Inverting 5yr land_ocean_edited with 1.2M verif pts at start, 230km gr supp3, 0,1,m2,k1,o1,a, zone 2050km, area_type 2, upcontver with disable_support 0"

//Split input file into pts that are used in inversion, and pts that aren't.
//0 - Do not split input file.  Use all input points in inversion.
//1 - Last "chunk" points are not used in inversion, only to verify results.
//2 - First "chunk" pts after "buffer" aren't used in inversion, only to verify results.
input.options.split_input = 0;
input.options.buffer= 0;
//input.options.chunk = 600000;//About 10% of 5-yr file.
//input.options.chunk = 1200000;//About 20% of 5-yr file.
input.options.chunk = 240000;//About 20% of 1/5 of 5-yr file.
//input.options.chunk = 3000000;//About 50% of 5-yr file, which had 6098802 pts.
//input.options.chunk = 3996631;//50% of 7-yr file (7751341 pts) but with 1 week buffer extra.
//input.options.chunk = 3170361;//50% of 5-yr file (6098802 pts) but with 1 week buffer extra.
//input.options.chunk = 4879042;//80% to match 1 of 5 pts decimation.
//input.options.chunk = 5488922;//90% to match 1 of 10 pts decimation.
//input.options.chunk = 5788922;//a little more for quicker tests.
//input.options.chunk = 75000;//About 50% of short2 input file.
//input.options.chunk = 3009158;//Makes 7 yr input have same # pts as 5 yr.
//input.options.chunk = 1351895;//Makes BAD 7 yr input-orbit have same # pts as 5 yr.
//input.options.chunk = 2551895;//Makes BAD 7 yr input-orbit have same # pts as 5 yr-20%.
//input.options.chunk = 1652539;//Makes "fixed" 7yr input-orbit have same # pts as 5 yr
//input.options.chunk = 2852539;//Makes "fixed" 7yr input-orbit have same # pts as 5 yr-20%
//input.options.chunk = 500;//About 50% of Hudson input file.

//Type of grid to be loaded.
//00  -   Custom grid, never the same twice so it always writes to an auto-rm'd ls cache folder.
//01  -   Equal-GRID, azimuthally symmetric, with cap over north pole.
//02  -   Equal-AREA, azimuthally symmetric, with cap over north pole.
//        This grid has no "gaps" in it which is accomplished by making the
//        halfwidths slightly smaller than designated in lat/lon_spacing.
//        As a result, the spacings entered here are going to be slightly
//        reduced in the actual grid.
//03  -   Based on 2, but gr_op.specified_lat is used to rotate around
//        the x axis. spec_lon is set to 0, most other types set spec_lat=0
//        !?IF TYPE 2 IS CHANGED TO NEW TOPO2 FUNCTION, NEED TO CHANGE THIS??
//        !?MUCH LATER: YES!! I JUST TRIED TO USE TOPO2 BUT IT CRASHES!
//10  -   Based on 2, but uses new topo routine to calc (inside bin ONLY)
//        wcover/mean_depth, sets larger dists to 1/2 circumference.
//11  -   Based on 10, shifted so specified_lat/lon has a pt directly on it.
grid.options.type = 10;

//Grid spacings are recorded in km, but equal-grid is usually specified in degrees.
grid.options.lat_spacing = 800;
//grid.options.lat_spacing = 1*deg2km;
grid.options.lon_spacing = grid.options.lat_spacing;

//Grid types 3,11, region type 3, analysis_list_choice 1 all use this pt.
//Lon is in 0..360 format. (WEST OR EAST?)
//grid.options.specified_lat = 48.0;//Can't rotate much more than this yet!
//grid.options.specified_lat = 59.54833333;//Yakutat, AK
//grid.options.specified_lon = 220.26500000;
//grid.options.specified_lat = 60.155;//Lerwick
//grid.options.specified_lon = 358.86;
//grid.options.specified_lat = 59.0;//Hudson Strait test boundary.
//grid.options.specified_lon = 292.0;
grid.options.specified_lat = 0.0;

//Used in grid type 0. Size of custom_lats needs to equal size of custom_lons!
//NOT TO BE CONFUSED WITH SPECIFIED_LAT/LON ABOVE! THESE AREN'T WRITTEN INTO
//THE GRID CACHE FILENAME, AND IN FACT CAN'T BE USED WITH ANY CACHE AT ALL!
double custom_lats[] = { 90.0 };
double custom_lons[] = { 180.0 };

//Type of area calculation for each mascon:
//0 - OBSOLETE Uses lat/lon_ halfwidth for approx. rectangular area.
//1 - OBSOLETE Uses analysis_options.disk_radius to calc circular areas. BAD IDEA!!!
//2 - Calcs area accounting for curvature of the spherical earth, assumes
//    grid is composed of at least one pole cap and all other
//    grid points are layed out in rows, with each grid point
//    at a particular lat having the exact same area.
//3 - Defines area of FES pts at each FES lat using area_type 2. Then every
//    FES pt is matched to exactly 1 mascon, and all the FES areas in each
//    mascon are summed. Intended for use when upward-continuing FES fields
//    such as truncation errors.
grid.options.area_type = 2;

//Determines which (if any) model to load.  This option is over-written
//if the analysis type does not call for a model.
//BE CAREFUL! THIS ISN'T NEEDED FOR COMPARE_INSITU...(), IT'S ONLY NEEDED
//FOR ANALYSIS TYPES THAT REQUIRE THE MODEL, LIKE 1001,1002,1003, ETC.
//IN FACT, IT SEEMS TO *SCREW UP* THE COMPARE_INSITU...() FUNCTION IF USED!!!
//0   - Do not load any model.
//1   - FES2004 amps/phases at gr pts using land pts too, and is_pt_in_bin.
//2   - Copy of 1, but using spherical_distance to resolve conflicts.
//3   - Copy of 2, but weighted using FES pts' areas.
//10  - Padman's AOTIM-5 tide model amps/phases using gaussian interp.
///////Compound model types - load 2 or more models and compare them.///////
//101 - FES2004 (type 1) minus AOTIM-5 (type 2) amps/phases
model_options.type = 0;

//Region code - divide grid up into local regions that will be analyzed
//separately.
//region_type can be:
//00   - Disables region code altogether.
//01   - N. Hemisphere equal area regions.
//02   - Entire grid is covered by one huge region.  Should work like 0, but
//       uses region code machinery so that create_least_squares() works.
//       Type 2 over-rides cap size, overlap, spacings, boundaries to make
//       sure that any least squares cache subfolder with type 2 is the same
//       as any other as long as they use identical grids.
//03   - Same as 1, shifted so reg.op.specified_lat/lon has a pt on it.
region_options.type = 2;

//For equal-grid region types, spacings are in units of degrees.
//For equal-area region types, spacings are in kilometers.
region_options.lat_spacing = 2700.0;
//region_options.lat_spacing = 3.0*grid.options.lat_spacing;//Useful for type 11.
//region_options.lon_spacing = region_options.lat_spacing;//Usually a good idea.
region_options.lon_spacing = 7000.0;//Sometimes regions need more longitude room.

//Used to determine number of edge points.
//region_options.overlap_distance = 0.0;
region_options.overlap_distance = 2850.0;
//region_options.overlap_distance=(double)((int)(3.0*grid.options.lat_spacing));

//This variable determines which of the defined regions will be analyzed.
//0 - Start at region 1 and work towards the very end.
//1 - Only analyze the reg which has (in center) reg_op.specified_lat/lon.
//    REGION TYPE 3 ALSO USES SPECIFIED_LAT/LON!!
//2 - Start at end and work backwards towards the beginning.
//3 - Start at end and work backwards, stop when region lat changes.
//    This creates output only near equator, where track density is low.
//    THIS OPTION NEEDS TO BE CHANGED WHEN S. GRID AND REGION ARE USED!
//4 - Analyze regions listed in manual_list (by number not index!) ONLY.
//5 - Analyze regions listed in manual_list first, then analyze all others.
region_options.analysis_list_choice = 0;
//Lon is in 0..360 format. (WEST OR EAST?)
region_options.specified_lat = grid.options.specified_lat;
region_options.specified_lon = grid.options.specified_lon;
//region_options.specified_lat = 80.0;
//region_options.specified_lon = 30.0;
long long manual_list[] = {1};//1st reg is 1, not 0.

//analysis_options.type controls the type of analysis performed.
//
//-1  - SKIP ANALYSIS
//00  - Load previously computed output using files referenced by
//      GAIAFOLDER/output/output_filenames.txt. Note: grid params are
//      over-ridden by grid parameters from the ones in output/options.txt
//
//INPUT ONLY - ACCELERATION SPACE
//01  - Each grid pt is the center of a bin, equal area or grid, with space
//      around it defined by lat(lon)_spacing.  Search the input file
//      for accels in each bin, and record the number of accels in the bin.
//      Also includes routines that are used to check region code.
//02  - Define bins as in type 1, examine input values inside each bin and
//      compute the RMS in each bin.
//03  - Each grid pt is the center of a bin, equal area or grid, with space
//      around it defined by lat(lon)_spacing.  Search the input file
//      for accels in each bin, and fit a constant term, secular trend,
//      and sinuisoids with periods from "parameter" using GSL QR decomp-
//      osition. Return the coefficients of the best-fit cosine and sine.
//04  - Same as 3, but returns amp and phase instead of cos and sin amps.
//05  - For use with region code only.  Make map where region "i" has the
//      value "i" in each of its center grid points.
//06  - For use with region code only.  Make map where each region has the
//      value "input.numpoints" in each of its center grid points.
//07  - For use w region code only.  Make 1 map per region where center pts
//      are 1.0 and ovlp pts are 0.5.
//8/9 - Evaluate phase of specified tidal constituents every time
//      GRACE flies directly over each bin (only 1st st pt in bin!)
//      Then output:
//      Type 8 -  average phase diff in degrees, rms about mean.
//      Type 9 -  # of phase diffs > phaseparam1, % of ".
//      THESE TYPES REQUIRE REGION CODE SO BININDICES IS AVAILABLE!
//10  - Evaluate phase of specified tidal constituents every time
//      GRACE flies directly over each bin (only 1st st pt in bin!)
//      Note that the ABSOLUTE phase is being recorded, not the
//      phase difference!
//      Then create a 2D plot for each bin showing how many times
//      each phase was found to have a certain value. The phases
//      are modulo 360 and grouped into chunks of size "lambda".
//      THIS TYPE REQUIRES REGION CODE SO BININDICES IS AVAILABLE!
//11  - Evaluate RELATIVE phase of specified tidal constituents every time
//      GRACE flies directly over each bin (only 1st st pt in bin!)
//      Note that the RELATIVE phase difference is being recorded, not the
//      absolute phase as in type 10!
//      Then create a 2D plot for each bin showing how many times
//      each phase difference was found to have a certain value. The phases
//      are grouped into chunks of size "lambda" up to 360.
//      THIS TYPE REQUIRES REGION CODE SO BININDICES IS AVAILABLE!
//12  - Copy of 10, but only makes plots for poles and one at each latitude.
//13  - Copy of 11, but only makes plots for poles and one at each latitude.
//14  - Maps first "regparamint1" input pts' positions as separate points.
//15  - Just like 11/13, but records time instead of phase.
//16  - Maps first "regparamint1" swath accels, with color scale.

//INPUT ONLY - MASCON SEA SURFACE HEIGHT SPACE
//102 - Invert input data via QR decomposition from gsl (~ least squares).
//      Returns coefficients of best-fit cosine and sine.
//103 - Same as 102, but returns amplitude and phase instead of cos and sin.
//104 - Invert input data via LU decomposition from gsl (~ least squares).
//      Returns coefficients of best-fit cosine and sine.
//105 - Same as 104, but returns amplitude and phase instead of cos and sin.
//106 - Invert input data via SVD from gsl.
//      Returns coefficients of best-fit cosine and sine.
//107 - Same as 106, but returns amplitude and phase instead of cos and sin.
//108 - Invert input data via SVD from gsl AFTER conversion to square Ax=b.
//      Returns coefficients of best-fit cosine and sine.
//109 - Same as 108, but returns amplitude and phase instead of cos and sin.
//110 - Invert input via truncated square SVD using jacobi (slow,accurat)
//      Returns coefficients of best-fit cosine and sine.
//111 - Same as 110, but returns amplitude and phase instead of cos and sin.
//112 - Invert input data via trunc. square SVD using mod (fast,!accurate)
//      Returns coefficients of best-fit cosine and sine.
//113 - Same as 112, but returns amplitude and phase instead of cos and sin.
//114 - Invert input data via sliced square SVD using jacobi (slow,accurat)
//      Returns coefficients of best-fit cosine and sine.
//      Slicing means SVD matrix is set to all zero except one sing value
//      that's chosen in the code in definitions.hpp.
//115 - Same as 114, but returns amplitude and phase instead of cos and sin.
//116 - Invert input data via LU decomposition after adding gaussian
//      smoothing term to the penalty function, lambda and regparam1 below!
//      Returns coefficients of best-fit cosine and sine.
//117 - Same as 116, but returns amplitude and phase instead of cos and sin.
//118 - Invert input data via LU decomposition after adding Tikhonov
//      regular. term to the penalty function, lambda needs a value below!
//      Returns coefficients of best-fit cosine and sine.
//119 - Same as 118, but returns amplitude and phase instead of cos and sin.
//120 - Invert input data via *REFINED* LU decomp after adding Tikhonov
//      regular. term to the penalty function, lambda needs a value below!
//      Returns coefficients of best-fit cosine and sine.  Uses convergence.
//121 - Same as 120, but returns amplitude and phase instead of cos and sin.
//122 - Invert input data via LU decomposition after INPLACE adding Tikhonov
//      regular. term to the penalty function, lambda needs a value below!
//      Returns coefficients of best-fit cosine and sine. Similar to 118.
//123 - Same as 122, but returns amplitude and phase instead of cos and sin.
//124 - Invert input data via Householder transformations after INPLACE adding Tikhonov
//      regular. term to the penalty function, lambda needs a value below!
//      Returns coefficients of best-fit cosine and sine. Similar to 122.
//125 - Same as 124, but returns amplitude and phase instead of cos and sin.
//
//MODEL ONLY.
//1001 - Map ALL model amps for model type listed in model_options.type.
//1002 - Map ALL model phases for model type listed in model_options.type.
//1003 - Create amp/phase output for model_options.type that can be loaded.
//
//MODEL AND input TOGETHER.
//2001 - input amplitude error (from type 4) divided by model ampl. (from
//       model type) NOT FINISHED - OR EVEN BEGUN, FOR THAT MATTER!
//
//SYNTHETIC - OUTPUT FILES USED FOR CREATING SIMULATED ACCELS, ETC.
//10001- Parameters - const, lin, m2, annual.  Amps on order of 1cm.
//       Const. - positive from 55-60, 65-70, 75-80,85-90, else negative.
//       Trend  - negative from ", positive otherwise. - only on land!
//       M2     - lat-dependent amplitude (highest at pole), phase = lon-dep
//                (constant in each quadrant) - only in ocean!
//       Annual - const amp in ocean only- phase varies smoothly with lon.
//       Saves cosine and sine coefficients.
//10002- Same as 10001, but saves amp and phase instead of cos and sin.
//10003- Parameters - const, lin, m2, annual.  Amps on order of 100cm.
//       Const. - positive from 50-70, negative otherwise
//       Trend  - negative from 55-60, 65-70, 75-80,85-90, positive
//                otherwise. - only on land!
//       M2     - lat-dependent amplitude (highest at pole), phase = lon-dep
//                (constant in each quadrant) - only in ocean!
//       Annual - const amp in ocean only- phase varies smoothly with lon.
//       Saves cosine and sine coefficients.
//10004- Same as 10003, but saves amp and phase instead of cos and sin.
//10005- Parameters - const, lin, m2, annual.  Amps on order of 10cm.
//       Const. - positive south of 75N, negative otherwise
//       Trend  - lon-dep: negative from 0-60, 120-180, 240-300, positive
//                otherwise. - only on land!
//       M2     - lat-dependent amplitude (highest at pole), phase = lon-dep
//                (constant in each quadrant) - only in ocean!
//       Annual - const amp in ocean only- phase varies smoothly with lon.
//       Saves cosine and sine coefficients.
//10006- Same as 10005, but saves amp and phase instead of cos and sin.
//10007- Parameters - const, lin, m2, annual.  Amps on order of 10cm.
//       Const. - single pt at 60N, 0E.
//       Trend  - single pt at 60N, 90E.
//       M2     - single pt at 60N, 180E.
//       Annual - single pt at 60N, 270E.
//       Saves cosine and sine coefficients.
//10008- Same as 10007, but saves amp and phase instead of cos and sin.
//10009- Parameters - const, lin, m2, annual.  Amps on order of 10cm.
//       All at single pt at 62N, 90E.
//       Saves cosine and sine coefficients.
//10010- Same as 10009, but saves amp and phase instead of cos and sin.
//10011- Parameters - const, lin, m2, annual.  Amps on order of 10cm.
//       All at single pt at NP.
//       Saves cosine and sine coefficients.
//10012- Same as 10011, but saves amp and phase instead of cos and sin.
//10013- Parameter - const only, 10cm in ocean only.
//       Saves cosine and sine coefficients.
//10014- Copy of 10012, but annual is replaced with K1.
//10015- Parameters - const, lin, m2, k1.  Amps on order of 10cm.
//sim15  All at one pt at 62N,35E - cent of rg #2, 3k rg, 1k cap,1250 ovlp
//sim61  Saves cosine and sine coefficients.
//sim73
//10016- Same as 10015, but saves amp and phase instead of cos and sin.
//10017- Copy of 10006, but annual is replaced with K1.
//10018- Copy of 10017, but K1 and trend patterns are swapped.
//       Parameters - const, lin, m2, k1. Amps on order of 10cm.
//       Const. - positive south of 75N, negative otherwise
//       Trend  - uniform trend in ocean only.
//       M2     - lat-dependent amplitude (highest at pole), phase = lon-dep
//                (constant in each quadrant) - only in ocean!
//       K1     - uniform amp on land only w lon-dep phase: -90 from 0-60,
//                120-180, 240-300, +90 otherwise.
//10019- Copy of 10018, but the constant is uniformly positive everywhere.
//10020- Constant has 10cm signal everywhere.
//10021- Copy of trend from 10017,10006.  2cm in amp.
//       Trend  - neg from lon=0-60, 120-180, 240-300, else pos. land only!
//10022- Copy of 10018, but const looks like a bullseye.
//sim25  Parameters - const, lin, m2, k1. Amps on order of 10cm.
//       Const. - positive between 70N and 80N, negative otherwise
//       Trend  - uniform trend in ocean only.
//       M2     - lat-dependent amplitude (highest at pole), phase = lon-dep
//                (constant in each quadrant) - only in ocean!
//       K1     - uniform amp on land only w lon-dep phase: -90 from 0-60,
//                120-180, 240-300, +90 otherwise.
//10023- Copy of 10022, but no trend.
//sim26  Parameters - const, m2, k1. Amps on order of 10cm.
//       Const. - positive between 70N and 80N, negative otherwise
//       M2     - lat-dependent amplitude (highest at pole), phase = lon-dep
//                (constant in each quadrant) - only in ocean!
//       K1     - uniform amp on land only w lon-dep phase: -90 from 0-60,
//                120-180, 240-300, +90 otherwise.
//10024- Copy of 10022, but M2 and K1 are swapped.
//sim27  Parameters - const, lin, m2, k1. Amps on order of 10cm.
//       Const. - positive between 70N and 80N, negative otherwise
//       Trend  - uniform trend in ocean only.
//       K1     - lat-dependent amplitude (highest at pole), phase = lon-dep
//                (constant in each quadrant) - only in ocean!
//       M2     - uniform amp on land only w lon-dep phase: -90 from 0-60,
//                analysis_options.type == 0120-180, 240-300, +90 otherwise.
//10025- Copy of 10023, but no M2 or const. K1 on land only.
//sim28  Parameters - k1. Amplitude is 2cm.
//       K1     - uniform amp on land only w lon-dep phase: -90 from 0-60,
//                120-180, 240-300, +90 otherwise.
//10026- Copy of 10022, but only trend and K1.
//sim29  Parameters - lin, k1. Amps on order of 10cm.
//       Trend  - uniform trend in ocean only.
//       K1     - uniform amp on land only w lon-dep phase: -90 from 0-60,
//                120-180, 240-300, +90 otherwise.
//10027- Copy of 10026, but only trend and K1 north of 65N.
//sim30  Parameters - lin, k1. Amps on order of 10cm.
//       Trend  - uniform trend in ocean north of 65N only.
//       K1     - uniform amp on land only w lon-dep phase: -90 from 0-60,
//                120-180, 240-300, +90 otherwise. (only north of 65N)
//10028- Copy of 10027, but trend replaced with const with same pattern.
//sim31  Parameters - const, k1. Amps on order of 10cm.
//       Const  - uniform trend in ocean north of 65N only.
//       K1     - uniform amp on land only w lon-dep phase: -90 from 0-60,
//                120-180, 240-300, +90 otherwise. (only north of 65N)
//10029- Copy of 10027, but only trend north of 65N.
//sim32  Parameters - lin. Amp = 10cm.
//       Trend  - uniform trend in ocean north of 65N only.
//10030- Copy of 10029, but constant instead of trend.
//sim33  Parameters - constant. Amp = 10cm.
//       Const - uniform value in ocean north of 65N only.
//10031- Copy of 10029, but rotated by 90 degrees towards the west.
//sim34  Parameters - lin. Amp = 10cm.
//       Trend  - uniform trend in ocean north of 65N only.
//
//10032- Secular trend in ocean as defined by wcover, uniform amplitude of 10cm.
//sim35
//
//10033- Constant in ocean, uniform amplitude of 10cm.
//sim36  - sim 36 is 230km, sim 53 is 500km eq area cap.
//
//10034- K1 Sine in ocean, uniform amplitude of 10cm.
//sim37
//
//10035- K2 Sine in ocean, uniform amplitude of 10cm.
//sim38
//
//10036- K1 Alias (RL2006) Sine in ocean, uniform amplitude of 10cm.
//sim39
//
//10037- K2 Alias (RL2006) Sine in ocean, uniform amplitude of 10cm.
//sim40
//
//10038- K2 Alias (RL2006) Sine in ocean north of 65N, uniform amp of 10cm.
//sim41
//
//10039- Const is +10cm between 70N and 80N, -10cm otherwise
//sim42
//
//10040- Const is +10cm between 55N and 65N, -10cm otherwise
//sim43
//
//10041- K1 +10cm between 55N and 65N, 0 otherwise.
//sim44
//
//10042- K2 +10cm between 55N and 65N, 0 otherwise.
//sim45
//
//10043- P1 +10cm between 55N and 65N, 0 otherwise.
//sim46
//
//10044- M2 +10cm between 55N and 65N, 0 otherwise.
//sim47, also sim68 with new grid.
//
//10045- K1 is +10cm between 55N and 65N (and 0E and 180E), 0 otherwise
//sim48
//
//10046- K1 is +10cm everywhere
//sim49
//
//10047- M2 is +10cm everywhere
//sim50
//
//10048- Constant is +10cm everywhere
//sim51
//
//10049- Trend  is +10cm/year everywhere
//sim52
//
//10050- Copy of 10023, but no constant, M2 and K1 are identical, sin/cos.
//sim62  Parameters - m2, k1. Amps = 10cm.
//       M2,K1  - phase = lon-dep (constant in each quadrant) - only in ocean!
//       Here ocean is defined using wcover > 0.
//10051- Copy of 10023, but no constant, M2 and K1 are identical, sin/cos.
//sim63  Parameters - m2, k1. Amps = 10cm.
//       M2,K1  - uniform amp on land only w lon-dep phase: -90 from 0-60,
//                120-180, 240-300, +90 otherwise.
//       Here ocean is defined using mascon elevation < 0.
//10052- Copy of 10050, but added trend along with identical M2, K1. sin/cos.
//sim64  Parameters - trend, m2, k1. Amps = 10cm.
//       Trend  - uniform trend on land.
//       M2,K1  - phase = (constant in quadrants NOT centered on NP - only in ocean!
//       Here ocean is defined using wcover > 0.
//10053- Copy of 10051, but added trend along with identical M2, K1. sin/cos.
//sim65  Parameters - trend, m2, k1. Amps = 10cm.
//       Trend  - uniform trend on land.
//       M2,K1  - uniform amp on land only w lon-dep phase: -90 from 0-60,
//                120-180, 240-300, +90 otherwise.
//       Here ocean is defined using wcover > 0.
//10054- Parameters - const, lin, m2, k1.  Amps = 10cm.
//sim74  Const. - single pt at 60N, 0E.
//       Trend  - single pt at 60N, 90E.
//       M2     - single pt at 60N, 180E.
//       K1     - single pt at 60N, 270E.
//       Saves cosine and sine coefficients.
//10055- Trend  is +10cm/year on land.
//       Here ocean is defined using wcover > 0.
//10056- Constant 10cm at one pt at 62N,35E.
//
//
//ANCILLARY PLOTS - ELEVATION, TEMPERATURE, ETC.
//20001- Elevation of grid point, in meters.  Negative values = ocean depth.
//20002- Pos elevations = 10, negative elevations = 0, to see coastlines better.
//
//INSITU PLOTS - LOCATIONS OF TIDE GAUGES, ETC.
//30001- Locations of tide gauges.
analysis_options.type = 20001;

//Sets regularization weight in many advanced types, beginning with 116.
//Also used in type 9 as a maximum phase (deg) difference between passes.
//Also used in type 10 as the size of each bin in the phase line plot.
analysis_options.lambda = 20;

//Sets a double float reg. parameter, set by the type chosen above.
//For instance, in 116,117, regparam1 is used as the gaussian halfwidth.
analysis_options.regparam1 = 78.0;
analysis_options.regparam2 = 1.0;//Not used at the moment.

//Sets a secondary integer reg. param, set by the type chosen above.
//Type 14 uses this to control how many input pts are plotted.
analysis_options.regparamint1 = 2;

//Refined inversions compare the initial misfit to the refined misfit and
//measure rms diffs.  Loop the refine code until rms diff / rms < conv.
analysis_options.convergence = 1E-15;

//Parameters (periods for sines, constant, trend, etc) used in analysis.
//               --See full list of defined periods in definitions.hpp.
//New parameter types should be added here.

//These arrays are coupled with matching array of parameter types:
//0    - Constant
//1    - Secular trend (Time^1)
//2    - Acceleration (Time^2)
//3-99 - Time^n polynomial term
//101  - Sinusoid (with period in seconds in parameter[]).

//GAIA can also solve for scale factors of arbitrary time series: 
//10001- Arbitrary secular trend (to test arbitrary inversions). 
//10101- Arbitrary cos(omega*t) term (to test arbitrary inversions). 
//10102- Arbitrary sin(omega*t) term (to test arbitrary inversions). 
//10201- S&P 500 daily closing values from 2002-2010 (to test arbitrary inversions). 
//11001- Daily Arctic OBP in cm water, highpass filtered, avg'd from Bering/NP gauges. 
//11002- ABPR1 hourly, de-trended, in cm water. 
//11003- ABPR1 hourly, de-trended, de-tided, in cm water. 
//11004- ABPR1 15 min, de-trended, in cm water. 
//11101- Upernavik (West Greenland) GRACE monthly (2002-2010) timeseries, detrended.
//12001- Use insitu time series (hopefully just one) as an arbitrary time series.

//Next type: 
//20001-  
//Note that the following p_types don't use the corresponding
//entries in the parameter array: 0-99,10000-20000(?)

//double parameter[] = {0,0,0};
//long long p_type[] = {0,1,11101};

//double parameter[] = {m2_period,k1_period,o1_period,p1_period};
//long long p_type[] = {101,101,101,101};

//double parameter[] = {mf_period,mf_period};
//long long p_type[] = {10101,10102};

//double parameter[] = {0,0};
//long long p_type[] = {11001,10001};

//double parameter[] = {0};
//long long p_type[] = {0};

double parameter[] = {1};
long long p_type[] = {1};

//double parameter[] = {p1_period};
//long long p_type[] = {101};

//double parameter[] = {m2_period};
//long long p_type[] = {101};

//double parameter[] = {1,m2_period};
//long long p_type[] = {1,101};

//double parameter[] = {k1_period};
//long long p_type[] = {101};

//double parameter[] = {m2_period,k1_period};
//long long p_type[] = {101,101};

//double parameter[] = {0,1,m2_period,k1_period};
//long long p_type[] = {0,1,101,101};

//double parameter[] = {1,m2_period,annual_period};
//long long p_type[] = {1,101,101};

//double parameter[] = {k1_period,k1_alias_period};
//long long p_type[] = {101,101};

//double parameter[] = {k2_period,k2_alias_period};
//long long p_type[] = {101,101};

//double parameter[] = {1,m2_period,annual_period};
//long long p_type[] = {1,101,101};

//double parameter[] = {0,1,2,3,11,m2_period};
//long long p_type[] = {0,1,2,3,11,101};

//double parameter[] = {0,1,m2_period,o1_period,k1_period,s2_period,p1_period};
//long long p_type[] = {0,1,101,101,101,101,101};

//double parameter[] = {0,1,2};
//long long p_type[] = {0,1,2};

//double parameter[] = {0,1,m2_period,k1_period,o1_period,annual_period};
//long long p_type[] = {0,1,101,101,101,101};

//double parameter[] = {0,1,m2_period,semiannual_period,annual_period};
//long long p_type[] = {0,1,101,101,101};

//double parameter[] = {m2_period,s2_period,k1_period,n2_period};
//long long p_type[] = {101,101,101,101};

//double parameter[] = {mm_period,semiannual_period,annual_period};
//long long p_type[] = {101,101,101};

//double parameter[] = {0,1,2,m2_period,k1_period,o1_period,p1_period,semiannual_period,annual_period};
//long long p_type[] = {0,1,2,101,101,101,101,101,101};

//double parameter[] = {0,1,q1_period,n2_period,mf_period,mm_period,semiannual_period,annual_period};
//long long p_type[] = {0,1,101,101,101,101,101,101};

//double parameter[] = {0,1,m2_period,k1_period,o1_period,p1_period,n2_period,semiannual_period,annual_period};
//long long p_type[] = {0,1,101,101,101,101,101,101,101};

//double parameter[] = {0,1,m2_period,k1_period,o1_period,annual_period};
//long long p_type[] = {0,1,101,101,101,101};

//double parameter[] = {0,1,m2_period,m4_period,k1_period,s2_period,o1_period,p1_period,n2_period,two_n2_period,k2_period,msqm_period,mtm_period,mf_period,mm_period,q1_period,semiannual_period,annual_period};
//long long p_type[] = {0,1,101,101,101,101,101,101,101,101,101,101,101,101,101,101,101,101};

//double parameter[] = {0,1,m2_period,k1_period,s2_period,o1_period,p1_period,n2_period,k2_period,q1_period,mf_period,mm_period};
//long long p_type[] = {0,1,101,101,101,101,101,101,101,101,101,101};

//double parameter[] = {m2_period,s2_period,k1_period,o1_period,s1_period,p1_period,n2_period,k2_period,mf_period,q1_period,m4_period,two_n2_period,mm_period,msqm_period,mtm_period};
//long long p_type[] = {101,101,101,101,101,101,101,101,101,101,101,101,101,101,101};

//double parameter[] = {0,1,m2_period,m2_plus_period,m2_minus_period,m4_period,k1_period,s2_period,o1_period,p1_period,n2_period,two_n2_period,k2_period,msqm_period,mtm_period,mf_period,mm_period,q1_period,lunar_sidereal_period,semiannual_period,annual_period};
//long long p_type[] = {0,1,101,101,101,101,101,101,101,101,101,101,101,101,101,101,101,101,101,101,101};

//The newer analysis types solve y=Hx+e, where "H" is set up as follows:
////Original, obsolete types:
//1   - GSL format, pt masses - obsolete, but contains kahan summation code.
//2   - GSL format, disk masses.  Uses nmax and disk_radius.
////Zone of influence types(uses zone value to choose ls cache subflder):
//101 - GSL format, pt masses.  Optimized for use with zone of influence.
////Supporting grid types (also uses gr support_type, mults for ls cache)
//201 - GSL format, pt masses.  Uses supporting grid and zone of influence.
//202 - GSL format, pt masses. Same as 201 but uses STL vectors in popultn.
analysis_options.h_matrix_choice = 202;

// 0 - Skip support grid creation.
// 1 - Create supporting grid using lat/lon_multiplier.
//     The multipliers determines how many support gr pts are associated
//     with each main gr pt. Uses support_lat/lon_multiplier values.
//     The supporting grid point code uses each main gr pt's lat/lon and
//     lat/lon halfwidth to break up each "quadrant" in each bin into
//     "lat_supp_mult*lonsupm" rectangular sections, then puts supp. gr pts
//     at corners that are not shared with other main grid points.
//     If either multiplier is set to 0, skip support grid creation.
//     If either mult. = 1, ERROR!, USED TO BE main gr pts have 4 supp pts: NW,SW,SE,NE.
//     If both multipliers are 2, all main gr pts have 8 supp points.
//     If both mults are 3, all main grid points have 24 supp gr pts, etc.
//     Each supp gr pt (and main gr pt) has a weight prop. to its area
//     so that closely-spaced support pts at northern edge of bin don't
//     weight the northern part of bin higher than southern edge which has
//     lower support grid point density.
grid.options.support_type=1;
grid.options.support_lat_multiplier=3;
grid.options.support_lon_multiplier=grid.options.support_lat_multiplier;

//Controls the vector over which a loop in create_least_squares()
//cycles through grid points. Thus, grid pts "too far" away from input pt
//don't have to be calculated in the contribution to the A,b matrices.
//0 - Go over all gr pts by creating trivial zone_of_influence vec of vecs.
//1 - Circularly symmetric zone with radius zone_radius (in km).
//NEXT ZONES ARE DEBUG, IDENTIFY NP AS I,J==0 (QUICK N DIRTY)
//2 - Same as 1, but excludes gr pt right below input pt (K1 investigation).
//3 - Same as 2, but NP st pts don't exclude NP mascon.
//4 - Same as 2, but ONLY NP st pts exclude NP mascon.
analysis_options.zone_of_influence_choice = 1;

//For zone_of_influence_choice 1, this is the farthest acceptable distance
//from input centerpoint to potential grid mascon before the grid's mass
//is assumed to affect the input's acceleration value.  In km.
//It's usually a good idea to set it equal to the region overlap distance.
//analysis_options.zone_radius = region_options.overlap_distance;
//analysis_options.zone_radius = 2850;//Standard for a long time.
analysis_options.zone_radius = 2050;//New standard.
//analysis_options.zone_radius = 1050;
//analysis_options.zone_radius = 100000;//Effectively infinite.

//Compare to in situ data (using the following definitions):
//autoFES   - Provides tide heights using automatic FES routines.
//TGs       - Tide gauge and bottom pressure recorder data.
//masconFES - Provides tide heights using mascon code, types 1003 or 0.
//real data - Mascon output interpolated to a location, making time series.

//0 - Don't run any comparison.
//1 - Compare autoFES to TGs.
//2 - Compare real data to "TGs - autoFES."
//3 - Compare masconFES to TGs.
//4 - Compare masconFES to autoFES, using TGs station locations and times.
insitu.comparison_type = 0;

//Verify how well mascons explain input accelerations by upward continuing
//mascons, then comparing original input file to resulting simulated input
//file.  input.options.split_input should probably be 1- that allows for
//input pts used in inversion to be compared separately from input pts not
//used in inversion.
//NOTE: This variable is copied to input.file_comparison_method in order to
//      control compare_input_files(input).
//0 - Do not perform upward continuation verification.
//1 - Verify by computing RMS of inv. pts (and, if present, verifc pts).
//    Same as 1 but non-tide soln performed separately, tides on land not used.
//    Also, tide/nontide noise analysis is performed for inver/verif accels.
//10- Just plot inversion and verification power spectra on the same plot.
analysis_options.upward_continuation_verification_choice = 0;
//If != 0, any verification choice calling for upward continuation will
//skip that process and instead load the subfolder "folder_for_skipping".
analysis_options.skip_upward_continuation_during_verification = 0;
analysis_options.folder_for_skipping = "/";

//Set to 1 if you don't want plots to be made in IDL.
plot_options.noplots = 0;

///////////////////////////////////////////////////////////////////////////
//Analysis options that are rarely- if ever- changed.
///////////////////////////////////////////////////////////////////////////

//0 - Does nothing.
//1 - Saves analysis results for each region, including overlap pts.
analysis_options.save_overlap_points=0;

//0 - Does nothing.
//1 - Saves all the input data found in each region.
analysis_options.save_regional_input_data=0;

//0 - Does nothing.
//1 - Immediately prints results at first grid point.
analysis_options.print_immediate_results=1;

//Pathnames of the stack input file (already processed by python script).
input.options.basefolder = "data/GRACE/";

//Topography data can be loaded from several sources:
//1  - Sean Swenson's topography file, with elevations and water depth.
//10 - BROKEN SRTM30+, global topo/bathy all in one big-endian binary file.
//11 - HUGE FILE SRTM30+, topo/bathy to 45N (5400 rows) all in one ascii file from Lin.
//20 - FES ascii file, with 999.9=land (except for NP). NO ELEVS/DEPTHS
ancillary.topo.type=1;

//Radius of cap over north pole, in km as of 2009-09-21.
//The cap size is usually specified by a flag for modern inversions, but
//old acceleration types like 1,2,3, or 4 reset this to 111 km.
//(Also, equal grid grids reset cap_size to 111 km as well.)
//"-1" is a flag that sets the cap_size to whatever size is necessary to
//match the area of the grid cap to match the areas of the nearest row of
//grid points. This requires adjusting the spacings in lat and lon,
//so they will change slightly.  UNLESS- if old accel plot is
//requested, autosizing is auto-deactivated and cap_size = 111 km.
grid.options.cap_size = -1.0;
//grid.options.cap_size = grid.options.lat_spacing;
//grid.options.cap_size = 2.0;
//cout<<"Grid auto cap sizing has been disabled for synthetic acceleration tests!!!!"<<endl;

//Specify grid type for final output (as opposed to internal grid which is always unstructured):
//0 - Unstructured grid output in main folder, no latlon grid output at all.
//1 - Unstructured grid output in main folder, latlon in "structured_output" subfolder.
//2 - Latlon in main folder, unstructured grid output in "unstructured_output" subfolder.
//3 - Latlon in main folder, no unstructured grid output at all.
//LEAVE ROOM FOR LATLON BASED ON SUPPORT GRID DENSITY!
grid.options.latlon = 0;
//When latlon grid's resolution is matched to original unstructured grid (for latlon > 0)
//the density of the latlon grid can by multiplied in both directions by these int factors (>=1)
grid.options.latlon_lat_multiplier=1;
grid.options.latlon_lon_multiplier=grid.options.latlon_lat_multiplier;

//If 1, analysis list created according to above choice is edited- the regs
//that have no lscache file present are not run.  If 0, run all listed regs.
//AMBIGUOUS: DO THEY NEED TO BE COMPLETED REGION CACHE FILES? FIND OUT!
region_options.analysis_list_cache_only = 0;

//Regional parallel processing (only used if region_options.type != 0):
//-1   - Autoselect.  Uses value of gaiafolder() to set in load_options().
//0    - No regional parallel processes- use sequential reg code if reg.np>0
//1    - Each reg. runs in a separate process, saves its results to a file.
//       These files are later combined into the final map.  The max num of
//       processes allowed to run simultaneously is "numprocesses".
region_options.parallel = 0;

//Temporal parallel processing:
//0    - No temporal parallel processing- run create_least_squares() once.
//1    - In create_least_squares_in_parallel(), split input into number of
//       chunks equal to numprocesses.  Populate A,b matrices using each
//       chunk in a different process using MPI.
analysis_options.time_parallel = 0;

//If ANY parallel code used, this option controls how many processes are run
//simultaneously.  This applies to regional parallelization and temporal
//paralleization (using MPI).
//-1   - Autoselect.  Uses value of gaiafolder() to set in load_options().
//                    This means it is reset every time options.txt is
//                    loaded, not just when options.txt is written.
//                    It looks for azonic, laptop strings and resets
//                    numprocesses to computer-specific values.
analysis_options.numprocesses = -1;

//Native models are usually hires compared to grid, so they're averaged.
//Controls Gaussian smoothed averages when loading model-for display 1001/2
//model_options.halfwidth=grid.options.lat_spacing/2;
//-1 is a signal to revert to the above value after grid is defined so that
//if a_o.type=0 and the loaded grid is different than specified here, it's
//based on the new grid rather than whatever's specified in this function.
model_options.halfwidth=-1;
//Each grid point averages the model points within this distance (in km).
//model_options.averaging_radius=2*model_options.halfwidth;
//Ditto.
model_options.averaging_radius=-1;

//For disk mascons, Legendre polynomials are expanded from n=0 to n=nmax.
analysis_options.nmax = 200;

//For disk mascons, each disk has this radius.  Also, OLD area_type routines
//used in define_grid described how many kg of water correspond to each
//cm in the output DID use this radius, but modern grids cover the earth, so
//disks aren't used for this purpose anymore because they leave gaps.
analysis_options.disk_radius = grid.options.lat_spacing/2.0;

//"boundaries" holds the limits of the grid, in terms of lat and lon.
//In general, the order is [min_lat, max_lat, min_lon, max_lon]
//Note: most applications requiring a local grid are best
//accomplished using a global grid (or at least a grid that covers
//the coordinates of the input data). Singling out a local area
//can be done using regions- simply choose a region size and
//specified_lat/lon to place it anywhere on the grid, then analyze
//only that region using whichever overlap is appropriate.
//The reason for this preference is that the input data is sorted
//according to which grid point it's closest to, so a grid that
//doesn't cover the entire area over which input data exists
//will simply pile all the "outside" input data into the outermost
//grid points. Using regions will avoid this problem.

//North pole to 49.8N.
grid.options.boundaries[0] =  49.8;
grid.options.boundaries[1] =  90.0;
grid.options.boundaries[2] =   0.0;
grid.options.boundaries[3] = 360.0;//*/

//Canadian TGs only.
/*grid.options.boundaries[0] =  57;
grid.options.boundaries[1] =  77;
grid.options.boundaries[2] =  360-130;
grid.options.boundaries[3] = 360-66;// */

//Greenland only.
/*grid.options.boundaries[0] = 58;
grid.options.boundaries[1] = 82;
grid.options.boundaries[2] = 285;
grid.options.boundaries[3] = 345;// */

//Around specified pt only.
/*double latbuffer = 1, lonbuffer = 1, jitter1=5, jitter2=0.1;
grid.options.boundaries[0] = grid.options.specified_lat-latbuffer-jitter1;
grid.options.boundaries[1] = grid.options.specified_lat+latbuffer+jitter2;
grid.options.boundaries[2] = grid.options.specified_lon-lonbuffer-jitter1;
grid.options.boundaries[3] = grid.options.specified_lon+lonbuffer+jitter2;// */

//Necessary to space regions using equal area code - radius of cap over
//north pole in km.
//"-1" is a flag that sets the cap_size to whatever size is necessary to
//match the area of the region cap to match the areas of the nearest row of
//region points. This requires adjusting the spacings in lat and lon,
//so they will change slightly.
//region_options.cap_size=deg2km;//1 deg. Nice for first region to be quick...
region_options.cap_size = 1700;
//region_options.cap_size = (grid.options.boundaries[1] - grid.options.boundaries[0])*deg2km;//Fill entire region with cap.
//region_options.cap_size = -1.0;

//Location of the tide gauge data.
insitu.tg.basefolder = "data/tide_gauges/";

//Location of the in situ cache folder.
insitu.cachefolder = "cache/insitu/";

//Location of grid cache folder.
grid.options.cachefolder = "cache/grid/";

//Location of model cache folder.
model_options.cachefolder = "cache/model/";

//Location of region cache folder.
region_options.cachefolder = "cache/region/";

//Location of input cache folder.
input.options.cachefolder = "cache/input/";

//Location of arbitrary cache folder.
input.arbitrary.cachefolder = "cache/arbitrary/";

//Location of least squares cache folder.
analysis_options.least_squares_cachefolder = "cache/least_squares/";

//Set to 1 to disable loading grid from cache.
//UPDATE: The grid computed from scratch (i.e. with no cache file or noload>0)
//        is slightly different than the grid loaded from cache.
//        Since 2010-08-26, ALL grids are stored to disk and reloaded, even
//        with noload>0. This insures that runs using a grid from scratch is
//        identical to a run using a grid loaded from a cache file.
//        Now all noload does is insure that the cache file is recomputed if it exists.
grid.options.noload = 0;

//Set to 1 to disable loading region from cache.
region_options.noload = 0;

//Set to 1 to disable loading model from cache.
model_options.noload = 1;//Does the model cache even WORK anymore?!

//Set to 1 to disable loading in situ data from cache.
insitu.noload = 1;//Needed because insitu cache doesn't work yet!

//Set to 1 to disable loading least squares matrices A,b from cache.
analysis_options.least_squares_noload = 0;

////////////////////////////////////////////////////////////////////////////
//Check variable length arrays for inconsistencies.
////////////////////////////////////////////////////////////////////////////

//Number of parameters needs to equal number of p_types.
long long num_parameter_types=(long long)(sizeof(p_type)/sizeof(*p_type));
long long num_parameters=(long long)(sizeof(parameter)/sizeof(*parameter));
if(num_parameters != num_parameter_types){
  cout<<"WARNING!! num_parameters = "<<num_parameters<<" while num_parameter_types = "<<num_parameter_types<<endl;
  grid.lat.resize(1); grid.lat.at(1) = 0;//Crash!
}

long long num_manual_list=(long long)(sizeof(manual_list)/sizeof(*manual_list));

//Number of custom grid lats needs to equal number of custom grid lons.
long long num_custom_lats=(long long)(sizeof(custom_lats)/sizeof(*custom_lats));
long long num_custom_lons=(long long)(sizeof(custom_lons)/sizeof(*custom_lons));
if(num_custom_lats != num_custom_lons){
  cout<<"WARNING!! num_custom_lats = "<<num_custom_lats<<" while num_custom_lons = "<<num_custom_lons<<endl;
  grid.lat.resize(1); grid.lat.at(1) = 0;//Crash!
}

////////////////////////////////////////////////////////////////////////////
//Create output folder, create scripts, write options to disk.
////////////////////////////////////////////////////////////////////////////

//The output folder is created here, and some scripts are saved to disk.
create_output_folder(argc,argv,plot_options);

write_options(options_file, input, region_options, ancillary, grid, model_options, analysis_options, plot_options, insitu, parameter, p_type, num_parameters, manual_list, num_manual_list, custom_lats, custom_lons, num_custom_lats);

FILE *output;

//Now that options file exists, copy to outputfolder.
string temp_string2 = "cp ";
temp_string2.append(options_file.c_str());
temp_string2.append(" \"");
temp_string2.append(plot_options.outputfolder.c_str());
temp_string2.append("\"");
//cout<<"options.txt copy string: "<<temp_string2<<endl;
output = popen(temp_string2.c_str(),"r");
pclose(output);

//Then copy source code archive (made by cgaia.sh) to the output folder.
temp_string2 = "cp ./c_source_files.zip \"";
temp_string2.append(plot_options.outputfolder.c_str());
temp_string2.append("\"");
//cout<<"zip source code string 2: "<<temp_string2<<endl;
output = popen(temp_string2.c_str(),"r");
pclose(output);
}

int main(int argc, char *argv[]){
  /**********************************************************************
  Purpose: This function is the grand overview of GAIA.
  ***********************************************************************
  Input:   argc        - number of words in command line string
           argv        - those words, used in output folder name
  ***********************************************************************
  Output:  Inverts data, verifies solution then writes results to disk.
  **********************************************************************/
  //Variable declarations- every variable type that ends in "_s"
  //is defined in GAIAFOLDER/c/definitions.hpp.
  long long i,k;//Counters.
  string options_file(gaiafolder());
  input_s input, unused_input;
  ancillary_s ancillary;
  grid_s grid, region, empty_region;
  model_s model;
  results_s results, regional_results,overlap_results;
  plot_options_s plot_options,temp_options,overlap_plot_options;
  insitu_s insitu;//For tide gauge data, etc.
  //These variables are used to display run time of program at end.
  time_t start_time,end_time;

  char s[max_length];
  string base;//This string forms the base part of the output filename.
  timespan_conversion_s timespan_conversion;//From secs to human units

  string temp_string;

  //Used to copy old options file instead of making it in create_options().
  string previous_options_file(gaiafolder());

  //Used to search the overlap folder for filenames of each region's files.
  DIR *cacheDIR;
  struct dirent *cacheDirEnt;
  string newfile;

  //////////////////////////////////////////////////////////////////////////
  //End of variable declarations, start of loading.
  //////////////////////////////////////////////////////////////////////////

  start_time = time(NULL);//Record time at start of program.

  //Location of the options file that create_options() creates
  //and load_options() loads, or the target of the copy operation that
  //loads a previously-generated options file.
  options_file.append("c/options.txt");

  //Options are chosen and written to disk in options_file.
  //argc and argv are used to create output folder.
  create_options(options_file, argc, argv);

  //Load options from file produced by create_options.txt
  load_options(options_file, input, region.options, ancillary, grid, model.options, results, plot_options, insitu,1);

  cout<<"Saving output to "<<plot_options.outputfolder<<endl;

  //Load ancillary data (topography, temperatures, etc) from disk.
  load_ancillary(ancillary, 1);

  //Next, a grid is defined for the inversion process in a "grid" struct.
  //Note: this is the "main" grid, which determines the spacing of points
  //used to create the map at the end.
  define_grid(grid, ancillary, results);

  //Now define the supporting grid.  The supporting grid is not used to
  //create the maps, and it does not affect the size of the A,B matrices.
  //Each main grid point has an (ostensibly arbitrary) number of supporting
  //grid points that are "tied" to it.
  //In upward continuation, the supporting grid pts are assigned the same
  //surface parameter value as its associated main grid point.
  //In Newtonian gravity inversions, supp gr pts partials are summed with the main
  //gr pt to form a single entry in the H matrix.
  define_support_grid(grid);

  //If desired, write supporting grid to disk so it can be plotted in IDL.
  if(1==2){
    write_supporting_grid(grid);
  }

  //Load in situ locations for comparisons.  Needs to be done before region
  //code so that optimized_region_indices code in define_region() can
  //prioritize region order based on which regions contain more stations.
  //(Although, frankly, I'm not sure I'll ever get around to coding this!)
  //model.options is passed because it contains averaging_radius and
  //halfwidth (used for Gaussian smoothing). These values are used to make
  //a list of nearby grid points and weighting factors for averaging the
  //residual signal at the tide gauges.
  if(insitu.comparison_type != 0) insitu_initialize(insitu, grid, model.options, ancillary);

  //If requested, define regions on which to solve for tides separately.
  region = define_region(region.options, grid, 1);

  //Now load JPL's input data - smoothed accelerations from GRACE A,B.
  if(input.options.noload == 0){
    load_input(input,grid,region);

    //Interpolate arbitrary time series to the input times.
    //If region code was used for loading, save the interpolated time
    //series to disk in cache/input.
    //This is done before split_input() so the interpolation happens
    //for verification accelerations as well.
    arbitrary_interpolate(input,region);

    //If requested, remove some input points from input and place them in
    //unused_input for later verification.  When regions are used as well,
    //unused_input is largely empty- the only change is that the original
    //input now has a new vector of uints called verification_points.  This
    //vector is the same length as binindices and is either 0 (pt used in
    //inversion) or 1 (pt not used in inversion, only for verification).
    if(input.options.split_input!=0) unused_input=split_input(input,region);
  }
  else{
    cout<<"Skipping input loading."<<endl;
  }

  //Now load the model if requested.
  load_model(model, grid);

  //Regions are not analyzed only in ascending order- 0,1,2,3. Instead,
  //they're ordered in different ways - region.options.analysis_list_choice
  //controls this. For instance, a single region can be analyzed on
  //its own. Needs to be called before least_squares_cache_initialize().
  if((long long)region.lat.size() != 0) build_analysis_region_list(grid, region);

  //If the analysis type is going to run create_least_squares(), then:
  //Search for matching least squares cache subfolder. If it doesn't exist,
  //create it and save name in results.options for future reference.
  //Also, copy options.txt to that folder so that a record is kept of which
  //parameters were used to create those matrices. NEED TO COPY GRID CACHE TOO,
  //IN CASE GRID IS DIFFERENT ON MACHINES THAT CREATE/USE LS CACHE FOLDER.
  if(results.options.uses_least_squares != 0 and (long long)region.lat.size()!=0)
  least_squares_cache_initialize(grid, region, input, results);
  //Otherwise set ls_cache_load to 0 so c_l_s() skips loading cache.
  else results.options.least_squares_cache_load = 0;

  //Since least_squares functions expect regions, warn if that's not true!
  if(results.options.uses_least_squares != 0 and (long long)region.lat.size()==0)
    cout<<"!!!WARNING!!! This analysis type ("<<results.options.type<<") calls for least squares, but region code is disabled because region.lat.size() == 0!"<<endl;

  //This forces ONLY center points to be copied. It's probably not needed,
  //but I'm making it explicit just to make sure that normal runs only copy
  //center points and not overlap points.
  results.options.copy_overlap_points=0;

  //Create an subfolder in the outputfolder, then make a new copy of
  //plot_options to deposit the new output files in the subfolder
  //rather than the outputfolder. It also stores the filenames of the
  //regional overlap files.
  if(results.options.save_overlap_points==1){
    //Create folder with desired foldername in the maps folder.
    plot_options.subfolder = plot_options.outputfolder;
    plot_options.subfolder.append("region_overlap_plots/");
    create_folder(plot_options.subfolder);

    overlap_plot_options = plot_options;
    overlap_plot_options.outputfolder = plot_options.subfolder;
    overlap_results.options = results.options;
    overlap_results.options.copy_overlap_points=1;
  }

  ///////////////////////////////////////////////////////////////////////
  //End of loading, start of analysis.
  ///////////////////////////////////////////////////////////////////////

  //If region.lat.size() == 0:
  //These commands analyze all input globally; they were primarily
  //used before I wrote the Newtownian gravity inversion routines.
  //The input data is analyzed using chosen params on the chosen
  //grid.  Additionally, a structure called "results.options" is
  //passed to perform_analysis() to specify the analysis method.
  //perform_analysis() modifies a structure "results" which
  //contains details of the analysis on the original grid. The
  //details may include things like fitted amplitudes of tidal lines, etc.
  //Also, type 7 needs regions defined but needs to be analyzed ONCE only.
  if((long long)region.lat.size() == 0 or results.options.type == 7){
    results.options.region_index = -1;//Note that reg code is off.
    results.options.prefix = "Entire grid";
    //Pass the real region struct bc it's needed for type 7.
    perform_analysis(input, grid, ancillary, model, results, plot_options, region, insitu);

    //If desired, print out 0th grid point results immediately.
    if(results.options.print_immediate_results==1) for(k=0;k<(long long)regional_results.outputs.size();k++) cout<<"Grid pt at lat "<<setiosflags(ios::fixed)<<setprecision(1)<<grid.lat.front()<<" and lon "<<setiosflags(ios::fixed)<<setprecision(1)<<grid.lon.front()<<" has "<<regional_results.titles.at(k)<<" = "<<setiosflags(ios::fixed)<<setprecision(4)<<regional_results.outputs.at(k).front()<<" "<<regional_results.units.at(k)<<endl;

  }
  //Else if region.lat.size() != 0:
  //Mult. regions will be analyzed separately using ONLY input data in each
  //geographical region, inverted to obtain solutions ONLY at grid points
  //in the region.  These solutions will be stored into regional_results
  //and "stitched" back together into a complete (global) solution in
  //results.
  //Region code can either be executed sequentially or in parallel.
  else{
    //Analyze regions either using sequential or parallel code:
    switch(region.options.parallel){
      default: cout<<"!!!!WARNING!!!!!! region.options.parallel "<<region.options.parallel<<" isn't recognized."<<endl;
      case 0:{//Sequential region code.
        sequential_regional_analysis(input, grid, ancillary, model, results, regional_results, overlap_results, plot_options, temp_options, overlap_plot_options, region, empty_region, insitu);
        break;
      }
      case 1:{//Parallel region code.
        parallel_regional_analysis(input, grid, ancillary, model, results, regional_results, overlap_results, plot_options, temp_options, overlap_plot_options, region, empty_region, insitu);
        break;
      }
    }//End of parallel processing switch-case-break.
  }//End of region.lat.size() != 0 else statement.

  //Calculate land tide before output is written to serve as an uncertainty estimate.
  estimate_land_tide(results,grid);

  //If requested, build latlon grid of "sufficient" density, copy output, then save to disk.
  copy_results_to_latlon(results, grid, plot_options);
  
  //Set the output filename base so that these output files are final.
  base = "plot_data";

  //Now that the analysis has been performed, it's time to output the data
  //in a format that the IDL program create_plots.pro can turn into a
  //plot.  Last variable is 1 because this is the final write_output.
  //Quick hack: if loading old output, don't output files to avoid filling
  //disk up w/ redundant files. Also don't output if analysis was skipped or latlon only wanted.
  if(results.options.type == 0 or results.options.type == -1 or grid.options.latlon == 3)
    cout<<"NOTE: not creating unstructured output files to avoid filling up disk with redundant files."<<endl;
  else
    write_output(base,results, grid, plot_options,1);

  //After outputting the final results, temp output files should be rm'd
  if((long long)region.lat.size()>0){
    for(i=0;i<(long long)temp_options.output_files.size();i++){
      if(remove(temp_options.output_files.at(i).c_str()) != 0) cout<<"Problem deleting temp file "<<temp_options.output_files.at(i)<<endl;
    }
  }//End of "if region.lat.size() > 0" statement.

  //If desired, save output from all points in each region,
  //center and overlap. This code creates the output_filenames.txt file.
  if(results.options.save_overlap_points==1){
    //Search the subfolder for files.
    cacheDIR = opendir(overlap_plot_options.subfolder.c_str());

    overlap_plot_options.just_the_filenames.clear();

    cacheDirEnt = readdir(cacheDIR);
    while(cacheDirEnt != NULL){
      strcpy(s,cacheDirEnt->d_name);
      newfile=s;
      //Only add files with 3 or more characters (i.e. not "." or "..").
      if((int)newfile.length() >= 3) overlap_plot_options.just_the_filenames.push_back(s);
      //cout<<"region overlap file: "<<s<<endl;
      cacheDirEnt = readdir(cacheDIR);
    }
    closedir(cacheDIR);

    //Sort the filenames alphabetically.
    sort(overlap_plot_options.just_the_filenames.begin(),overlap_plot_options.just_the_filenames.end());

    //Create output_filenames with those filenames.
    create_output_filenames(overlap_plot_options,results);

  }//End of "save_overlap_points".

  ///////////////////////////////////////////////////////////////////////
  //End of analysis, start of verification of results.
  ///////////////////////////////////////////////////////////////////////

  //Load in situ data for comparison with analysis results and/or model.
  if(insitu.comparison_type != 0) load_insitu(insitu,ancillary);

  //Code for FES 2D plots.
  if(1==2){
    //Plot accels vs time, as well as masconFES height values for the bin
    //where the accel was recorded, as well as interpolated masconFES height
    //values at the accel's ACTUAL location to test TG interpolation code, and
    //autoFES height value right underneath the satellite.
    //Plots first "regparamint1" swaths.
    //Needed by write_output().
    results_s temp_results;
    string xybase = "xyplotdata";
    grid_s temp_grid=grid,temp_region=region;
    plot_options_s temp_plot_options = plot_options;
    ancillary_s temp_ancillary;
    insitu_s temp_insitu,empty_insitu;
    input_s temp_input;
    //Used to push_back vectors of dbl vectors
    vector<double> temp_v_double;
    //push_back vs of vs of dbl vectors!
    vector< vector<double> > temp_vv_double;
    //push_back vs of string vectors!
    vector<string> temp_v_string;
    
    //Holds indices of pts in current "swath candidate"
    vector<long> current_swath_indices;

    //Override the usual map output.
    temp_results.options.output_choice = 2;

    //Because otherwise generate_input_timeseries doesn't use land mascons.
    for(i=0;i<(long long)temp_grid.wcover.size();i++) temp_grid.wcover.at(i) = 1.0;

    temp_input.options = input.options;
    temp_region.options.type = 2;
    temp_region = define_region(temp_region.options, temp_grid, 1);
    load_input(temp_input,temp_grid,temp_region);

    /////////////////////////////////////////////////
    //Input accel vs. masconFES of closest mascon.
    /////////////////////////////////////////////////
    
    //Loop until numswaths have been recorded or end of input is reached.
    i=0;//No swaths recorded yet for this plot.
    long long j=1;//Starting on second input entry so that c_s_i has at least 1 entry.
    current_swath_indices.push_back(0);
    while(i<results.options.regparamint1 and j<(long long)input.time.size()){
      //Is the next point too far away to be in the same swath, or is there NO next point?
      if(input.time.at(j) - input.time.at(current_swath_indices.back()) > 5 or j==(long long)input.time.size()-1){
        long numpoints = current_swath_indices.size();
        //If so, have enough points been found for this to be a swath?
        if(numpoints >= 240){
          
          //This is a brand new plot, so push_back using temp vectors/vvs.
          temp_results.xy.x_values.push_back(temp_vv_double);
          temp_results.xy.y_values.push_back(temp_vv_double);
          temp_results.xy.legends.push_back(temp_v_string);

          temp_results.xy.legends.back().push_back("Scaled Accels");
          temp_results.xy.titles.push_back(temp_results.xy.legends.back().back());
          temp_results.xy.titles.back().append(" vs ");
          temp_results.xy.legends.back().push_back("Binary MasconFES");
          temp_results.xy.titles.back().append(temp_results.xy.legends.back().back());
          sprintf(s,", swath #%3lld",i+1);
          temp_results.xy.titles.back().append(s);
          temp_results.xy.x_units.push_back("seconds");
          temp_results.xy.y_units.push_back("cm");

          int numlines=temp_results.xy.legends.back().size();
          temp_results.xy.x_values.back().resize(numlines);
          temp_results.xy.y_values.back().resize(numlines);

          for(k=0;k<numlines;k++){
            temp_results.xy.x_values.back().at(k).resize(numpoints);
            temp_results.xy.y_values.back().at(k).resize(numpoints);
          }
          
          //Save individual XY values for first line.
          for(k=0;k<numpoints;k++){
            temp_results.xy.x_values.back().at(0).at(k) = input.time.at(current_swath_indices.at(k))-input.time.at(current_swath_indices.front());
            temp_results.xy.y_values.back().at(0).at(k) = input.accel.at(current_swath_indices.at(k));
          }

          //Create fake insitu data from input coords, times.
          temp_insitu=empty_insitu;//To avoid problems caused by repeatedly running this code.
          temp_insitu.tg.times.resize(numpoints);
          temp_insitu.tg.sealevels.resize(numpoints);
          temp_insitu.tg.weighting_factors.resize(numpoints);
          temp_insitu.tg.nearby_grid_indices.resize(numpoints);
          temp_insitu.tg.nearby_grid_distances.resize(numpoints);
          for(k=0;k<numpoints;k++){
            temp_insitu.tg.titles.push_back("temp");
            temp_insitu.tg.grid.lat.push_back(input.lat.at(current_swath_indices.at(k)));
            temp_insitu.tg.grid.lon.push_back(input.lon.at(current_swath_indices.at(k)));
            temp_insitu.tg.times.at(k).push_back(input.time.at(current_swath_indices.at(k))+input.epoch);
            temp_insitu.tg.sealevels.at(k).push_back(1.0);
            temp_insitu.tg.weighting_factors.at(k).push_back(1.0);
            //cout<<"Made it.8 temp_input.binindices.size() = "<<temp_input.binindices.size()<<endl;
            temp_insitu.tg.nearby_grid_indices.at(k).push_back(temp_input.binindices.at(current_swath_indices.at(k)));
            temp_insitu.tg.nearby_grid_distances.at(k).push_back(1.0);
          }

          generate_input_timeseries(temp_insitu,results,temp_grid,0);

          //Save individual XY values for second line.
          for(k=0;k<numpoints;k++){
            temp_results.xy.x_values.back().at(1).at(k) = temp_results.xy.x_values.back().at(0).at(k);
            temp_results.xy.y_values.back().at(1).at(k) = temp_insitu.tg.input_sealevels.at(k).at(0);
          }
          i++;
        }
        //Either way, the last swath is gonna get cleared.
        current_swath_indices.clear();
        current_swath_indices.push_back(j);
      }
      else{//If index j is in the current swath, add it to the list.
        current_swath_indices.push_back(j);
      }
      j++;//Next input point.
    }//End of loop through input.
    scale_plot(temp_results,0,-1);//Scale line at first index by lines in second index (-1 means all others)
    write_output(xybase, temp_results, temp_grid, temp_plot_options, 1);
  }

  //CHECK THESE OUTPUT FILES TO MAKE SURE THAT FES IS BEING AVERAGED THE SAME USING 50KM AND 230KM
  //MASCON SPACING!
  if(1==2){
    long long j;
    results_s temp_results=results;
    temp_results.options.adjust_phases_choice = 1;
    adjust_phases(input,temp_results);
    convert_results_from_ampphase_to_sincos(temp_results);
    for(i=0;i<(long long)temp_results.outputs.size();i++){
      double kg_sum=0,pos=0,neg=0;
      for(j=0;j<(long long)temp_results.outputs.at(i).size();j++){
        kg_sum += temp_results.outputs.at(i).at(j)*grid.cm2kg.at(j);
        if(temp_results.outputs.at(i).at(j)>0) pos += temp_results.outputs.at(i).at(j)*grid.cm2kg.at(j);
        else neg -= temp_results.outputs.at(i).at(j)*grid.cm2kg.at(j);
        //cout<<setw(6)<<j<<", lat "<<setw(6)<<fixed<<setprecision(2)<<grid.lat.at(j)<<", lon "<<setw(6)<<fixed<<setprecision(2)<<grid.lon.at(j)<<", output: "<<setw(21)<<scientific<<setprecision(14)<<temp_results.outputs.at(i).at(j)<<", cm2kg: "<<setw(21)<<scientific<<setprecision(14)<<grid.cm2kg.at(j)<<endl;
      }
      cout<<"Parameter "<<temp_results.titles.at(i)<<" has kg_sum = "<<scientific<<setprecision(10)<<kg_sum<<endl;
      cout<<"Positive sum = "<<scientific<<setprecision(10)<<pos<<endl;
      cout<<"Negative sum = "<<scientific<<setprecision(10)<<neg<<endl;
    }
  }

  //Report the sum of the squares of tidal cos/sin coefficients over land and ocean
  //using wcover in order to divide the (accel) variance drop due to land/ocean tides 
  //by these values. Is the variance drop due to land tides smaller than that due to 
  //ocean tides even when normalized by the sum of the squares of the cos/sin coefs.
  //Don't divide by the number of ocean/land points, but weight by the area of each mascon.
  if(1==1){
    long long j;
    results_s temp_results=results;
    temp_results.options.adjust_phases_choice = 1;
    adjust_phases(input,temp_results);
    convert_results_from_ampphase_to_sincos(temp_results);
    double alltides_land_sum=0.0, alltides_ocean_sum=0.0;
    for(i=0;i<(long long)temp_results.outputs.size();i++){
      if(temp_results.options.output_tide.at(i) == 1){
        double land_sum=0.0,ocean_sum=0.0;
        for(j=0;j<(long long)temp_results.outputs.at(i).size();j++){
          if(grid.wcover.at(j) > 0.0) ocean_sum += pow(temp_results.outputs.at(i).at(j)*grid.cm2kg.at(j),2);
          else land_sum += pow(temp_results.outputs.at(i).at(j)*grid.cm2kg.at(j),2);
        }
        alltides_land_sum  += land_sum;
        alltides_ocean_sum += ocean_sum;
        cout<<"Parameter "<<temp_results.titles.at(i)<<" has land_sum "<<scientific<<setprecision(10)<<land_sum<<" and ocean_sum "<<scientific<<setprecision(10)<<ocean_sum<<endl;
      }
    }
    cout<<"Total land_sum = "<<scientific<<setprecision(10)<<alltides_land_sum<<" and total ocean_sum = "<<scientific<<setprecision(10)<<alltides_ocean_sum<<endl;
  }
  
  if(1==2){
    //Maps accelerations due to each mascon for the first "regparamint1"
    //accels, with color scale.
    results_s init_results, accels_results,temp_results=results;
    plot_options_s accels_plot_options, all_maps;
    input_s accels_input;
    string accels_base;

    //Resize arrays, define title and output folder and filename base.
    accels_results.titles.clear();
    accels_results.titles.push_back("Temp");
    accels_results.outputs.clear();//Only 1 output file at a time.
    accels_results.outputs.resize(1);//Only 1 output file at a time.
    accels_results.outputs.front().resize(grid.lat.size(),0.0);
    accels_results.options.parameter.clear();
    accels_results.options.parameter_type.clear();
    accels_results.options.parameter.push_back(0);
    accels_results.options.parameter_type.push_back(0);
    accels_results.options.sincos = 0;
    accels_results.options.uses_least_squares = 1;
    accels_results.options.output_type.clear();
    accels_results.options.output_type.resize(1);
    accels_results.options.output_tide.clear();
    accels_results.options.output_tide.resize(1);
    accels_plot_options = plot_options;
    all_maps = accels_plot_options;
    accels_results.options.output_choice=1;
    accels_results.options.current_phase_format = 0;
    init_results = accels_results;

    temp_results.options.adjust_phases_choice = 4;
    adjust_phases(input, temp_results);
    //Because otherwise perform_upward_continuation doesn't use land mascons.
    grid_s temp_grid=grid;
    for(i=0;i<(long long)temp_grid.wcover.size();i++) temp_grid.wcover.at(i) = 1.0;

    //Loop until numaccels have been recorded or end of input is reached.
    long long j=0;
    while(j<temp_results.options.regparamint1 and j<(long long)input.time.size()){
    
      accels_results = init_results;
      accels_input = input;
      //cout<<"Before upcont, accel["<<j<<"] = "<<scientific<<setprecision(3)<<accels_input.accel.at(j)<<endl;
      perform_upward_continuation2(accels_input, temp_grid, temp_results, accels_results,j);
      //perform_upward_continuation2(accels_input, temp_grid, temp_results, accels_results,-1);
      //perform_upward_continuation(accels_input, temp_grid, temp_results);
      //accels_input.options.output_filename = "input_and_orbit_gaia_global_region2.txt";
      //write_input(accels_input);
      //cout<<"After upcont, accel["<<j<<"] = "<<scientific<<setprecision(3)<<accels_input.accel.at(j)<<endl;

      //Calculate sum of all accelerations at this time, to compare to ORIGINAL input.accel.
      double accel_sum=0;
      for(k=0;k<(long long)accels_results.outputs.front().size();k++) accel_sum += accels_results.outputs.front().at(k);
      //cout<<"Accel_sum = "<<scientific<<setprecision(3)<<accel_sum<<endl;

      //Save misc info for this plot.
      sprintf(s,"Time: %lld, Sum: %.2e",accels_input.time.at(j)+accels_input.epoch,accel_sum);
      accels_results.titles.clear();
      accels_results.titles.push_back(s);
      sprintf(s,"mascon_accels_%03lld",j+1);
      accels_base=s;
      accels_results.units.push_back("m/s^2");
      
      write_output(accels_base, accels_results, temp_grid, accels_plot_options, 0);
      //Save the filenames for output_filenames (after loop).
      for(k=0;k<(long long)accels_plot_options.just_the_filenames.size();k++){
        all_maps.just_the_filenames.push_back(accels_plot_options.just_the_filenames.at(k));
      }
      j++;//Next input point.
    }//End of loop through input.
    //Create output_filenames with ALL those filenames.
    create_output_filenames(all_maps,accels_results);
  }

  //Load in situ data and compare it to the model or the input-derived
  //mascon signal. Note that model isn't sent, only model.options. That's
  //because the gaia option to load a model isn't used here- the model is
  //loaded from scratch using different functions.
  if(insitu.comparison_type != 0) insitu_comparison_verification(insitu, results, grid, model.options, plot_options);

  //Now verify output mascons by upward continuing them to form a simulated
  //input file and comparing it to the original input file.
  //
  //If results.options.upward_continuation_verification == 0, skip.
  //
  //If upward_continuation_verification > 0, then 1 (possibly 2) RMSs will
  //be computed- the RMS of input pts used for inversion (see split_input)
  //and the RMS of the input pts not used in the inversion.
  upward_continuation_verification(grid, input, unused_input, results, region, plot_options);

  end_time = time(NULL);//Record time at end of program.
  timespan_conversion.seconds = (long long)(end_time - start_time);
  cout<<"GAIA v"<<fixed<<setprecision(3)<<GAIA_VERSION<<" finished after"<<sec2human(timespan_conversion)<<"."<<endl;

  return 0;
}
